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ABSTRACT 

Parker's model is one of the most discussed mechanisms for coronal heating and has generated 
much debate. We have recently obtained new scaling results in a two-dimensional (2D) version of this 
problem suggesting that the heating rate becomes independent of resistivity in a statistical steady 
state [Ng and Bhattacharjee, Astrophys. J., 675, 899 (2008)]. Our numerical work has now been 
extended to 3D by means of large-scale numerical simulations. Random photospheric footpoint motion 
is applied for a time much longer than the correlation time of the motion to obtain converged average 
coronal heating rates. Simulations are done for different values of the Lundquist number to determine 
scaling. In the high-Lundquist number limit, the coronal heating rate obtained so far is consistent 
with a trend that is independent of the Lundquist number, as predicted by previous analysis as 
well as 2D simulations. In the same limit the average magnetic energy built up by the random 
footpoint motion tends to have a much weaker dependence on the Lundquist number than that in 
the 2D simulations, due to the formation of strong current layers and subsequent disruption when 
the equilibrium becomes unstable. We will present scaling analysis showing that when the dissipation 
time is comparable or larger than the correlation time of the random footpoint motion, the heating 
rate tends to become independent of Lundquist number, and that the magnetic energy production is 
also reduced significantly. 

Subject headings: magnetic reconnection — magnetohydrodynamics (MHD) — Sun: corona — Sun: 
magnetic topology 
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1. INTRODUCTION 

The enormous energy content of high-beta photo- 
spheric plasma flows has long been suggested as the 
source of energy that ultimately heats the million de- 
gree solar corona. Unambiguously identifying the ex- 
act mechanisms that transfer this kinetic energy into the 
overlying solar atmosphere and the exact nature of how 
the coronal magnetic field responds and converts this en- 
ergy into heat remains one of the long-standing issues in 
astrophysics. 

In this paper we investiga te an idealized model of the 
corona proposed by Parker ( 1972[ ) which applies to closed 
magnetic field structures whose field lines are embed- 
ded at both ends in the solar surface. Neglecting the 
curvature of the magnetic field, the corona is modeled 
in Cartesian geometry where an initially uniform mag- 
netic field along the Cz direction is line-tied at z = 
and z = L in perfectly conducting end-plates repre- 
senting the photosphere. Parker suggests that slow and 
continuous random shuffling of the footpoints at these 
end-plates, representing the turbulent buffeting of the 
coronal field embedded in the convecting photosphere, 
can tangle the field into a braided structure of sufficient 
complexity such that it cannot settle into a continuous 
smooth equilibrium, but rather necessarily evolves to one 
with tangential discontinuities. Whether or not true cur- 
rent singularities (as opposed to current layers with finite 
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thickness) can form in this scenario and whether or not 
continuous footpoint mappings necessarily imply a non- 
smooth topology has been the subject of intense debate 
in the decades that have passed since Parker's seminal 
proposal. Extended discussion of this matter is beyond 
the scope of the pr esent analysis, but it is appr opriate to 
reiterate here (c.f. Ng & Bhattacharjee 1998) that this 
question is not merely of academic interest. That the 
plasma gradients will tend towards singularities has im- 
portant bearing on the physics of magnetic reconnection 
and turbulence dyna mics in the corona. The interested 
reader is referred 



to 



3w| 



Ng fc B hattacharjee (1998), Low| 
( |2010| , [Huang et al.l ( |2010[ ), |Janse, Low,&: Parker, ( ,201 
and references therein. 

In a process Parker calls "topological dissipation", it 
is at these tangential discontinuities where the corona's 
small but ultimately finite resistivity induce the forma- 
tion of current sheets where magnetic energy is dissi- 
pated to heat the coronal plasma, and where magnetic 
reconnection proceeds to reduce the topological complex- 
ity of the coronal magnetic field. This esse ntial concept 



was further dev eloped in a series of studies ( Parker|"l979 
1983a|b 1988 |1994) and has become known colloqui 



ally as the "nanofiare model" of coronal heating. The 
appellation derives from the isolation of 10^^ erg flares 
as the constitutive energy release events which occur in 
"storms" of sufficient ferocity to heat the corona and 
adequately account for observed conductive and radia- 
tive losses. While the concept of topological dissipa- 
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tion can be seen as the prototypical "DC" mechanism 



for c oronal heating (c.f. Klimchuk 2006 Aschwanden 



2004 1 , the solar atmosphere surely admits more complex 



magnetic topology than is treated by the Parker model 
In fact, many investigators have pursued reconnection- 
based heating mechanisms using geometries that i nclude 



separators, separatrices and magnetic-nulls (see Priest 
et al.||2005| and references therein), and more recently 



by analyzing the magnetic topology of active regions ob 



served by Transition Region and Coronal Explorer (Lee 
et al |2010[ ). It remains clear however, that coronal loops 



are the basic building block of the solar corona, and that 
their examination first as isolated entities is crucial in 
laying the foundations for a b roader under standing of 
the corona and its activity (see Reale 2010 for a recent 
review) . 

The work we present her e is motivated by a recent 
study ( |Ng fc Bhattacharjee|2008^ hereafter NB08) which 
developed a simplified version of the Parker scenario. 
NB08 referred to this version as a constrained tectonics 
model 



at the 



( Priest et ah] 2002) because the random braiding 
ine-tied ends was restricted to depend on only 
one coordinate transverse to the initial magnetic field. 
This strong assumption has the consequence that it en- 
ables us to describe the complete dynamics of the sys- 
tem by a simple set of differential equations which are 
easily amenable to analytical and numerical solutions for 
prescribed footpoint motion. The geometric constraints 
imposed by our assumption preclude the occurrence of 
nonlinear effects such as reconnection and secondary in- 
stabilities, but enables us to follow for long times the dis- 
sipation of energy due to the effects of resistivity and vis- 
cosity. Using this model, it was shown both numerically 
and by scaling analysis that as long as the correlation 
time of turbulent photosphcric flow (tc) is much smaller 
than the characteristic resistive timescales {rf/), ohmic 
dissipation becomes independent of resistivity (?;). The 
absence of nonlinear effects in this model allows the per- 
pendicular magnetic field {B±) to grow to unphysically 
large values and is found to scale as 77"^/^. Furthermore, 
NB08 conjectured, by means of a heuristic scaling ar- 
gument, that even in the presence of reconnection and 
secondary instabilities, the heating rate would remain 
insensitive to resistivity. It is this conjecture that we ex- 
amine here using three-dimensional (3D) reduced MHD 
(RMHD) simulations. 

The Parker model has been s tudied exteiisi vely using 
3D MHD numerical s imulations (pVIikic et al.|1 989; Long- 
cope fc Sudan"1994Vii:inau di et al.|l l996', 'He ndrix et al.| 
1996; Galsgaard fc Nordlu nd||1996[ [Dm itruk et al.||1998[ 
Uomez et al.||2000[ |Rappazzo et al ||2008t [2010]among 
others.) Here we are interested in the precise scaling of 
dissipation with respect to plasma r esistivity. Our study 



is mos t similar in design to that of lLongcope fc Sudan 
( 1994) who used RMHD to simulate Parker's model with 



Lundquist numbers spanning an order or magnitude. In 
this range they found that both heating rate and perpen- 
dicular field production scales as ri~^^^. These numeri- 
cal results agreed with analysis based on Sweet-Parker 
reconnection theory and measurements of current sheet 
statistics. 

In this pape r, we recover the sca l ings o f heating rate 
and of B± of Longcope fc Sudan (1994) in the range 



they examined, and extend their results to a range of 
lower 77, where the results support a slower growth of B± 
which roughly scales as ri~^^^ and a heating rate that 
becomes insensitive to 77. We also demonstrate by simple 
scaling analysis that the transition between these scaling 
behaviors results from the diminishing effects of random 
photospheric motion as the energy dissipation timescale 
te becomes much smaller than the correlation time Tc, 
in accordance with NB08. 

The paper is organized as follows. In Section [2] we de- 
scribe the RMHD coronal heating model and numerical 
scheme. Section |3] gives the details of the simulation re- 
sults. Section |4] presents a scaling analysis describing the 
transition in the observed scaling behavior. Section |5] 
summarizes our conclusions to date and discusses out- 
standing physical issues. 

2. CORONAL HEATING MODEL AND NUMERICAL SETUP 

We assume that the coronal plasma is sufficiently low- 
beta so that the dynamics can be described by the 
RMHD equations. Many numerical studies of Parker's 
model mentioned above (wit h the notable exc eptions of 
Galsga ard fc Nordlund 1996; 'Gudiksen fc Nordlund 20051 
[Peter et al. 2006) arc based on the RMHD approxima- 
tion! ^I'hc RMHD equations are a simplified version of 
MHD applicable to systems where the plasma is domi- 
nated by a strong guide field such that the timescales of 
interest are slow compared with the characteristic Alfven 
timescale ta- These restrictions also imply incompress- 
ibility (V • v = 0) and the exclusion of magnetosonic 
modes (leaving only the shear Alfven modes propagat- 
ing in iz). The RMHD equations were first derived for 
the stu dy o f tokamak plas mas by |Kadomtsev fc Pogutse| 
(197A) and iStraussI (|1976|), which can be written in di- 
mensionless form as 



Ik 



8 T 



OA 



dz 



d(j) 
dz 



(1) 



(2) 



where A is the flux function so that the magnetic field 
is expressed as 6 = 62+ Bj^ = -I- V ^A x e^; (f) is 
the stream function so that the fiuid velocity field is ex- 
pressed as V = Vj^(?!)X Gz; Q. — is the z-component 
of the vorticity; J — —V\A is the z-component of the 
current density; and the bracketed terms are Poisson 
brackets such that, for example, [0, ^] = 4iyA,j. — (pxAy 
with subscripts here denoting partial derivatives. The 
normalized viscosity v is the inverse of the Reynolds num- 
ber Rv, and resistivity r] is the inverse of the Lundquist 
number S. The normalization adopted in Equations ([I]) 
and ([2]) is such that the magnetic field is in the unit of 
Bz (assumed to be a constant in RMHD); velocity is in 
the unit of va — Bz/{4:TTp)^^^ with a constant density 
p; length is in the unit of the transverse length scale 
Lj_; time t is in the unit of Lj_/va', rj is in the unit of 
AttvaL^l/c^] and v is in the unit of pvaL±. 

In this paper we i nves tiga te an idealized model of the 
corona proposed by IParkeri ((1972). In Parker's model, a 



solar coronal loop islreated as a straight ideal plasma 
column, bounded by two perfectly conducting end-plates 
representing the photosphere. The footpoints of the mag- 
netic field in the photosphere are frozen (line-tied). Ini- 
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tially, there is a uniform magnetic field along the z di- 
rection. The footpoints on the plates z = and z — L 
are subjected to slow, random motion ^(z — 0, t) and 
(j){z = L, t) that deform the magnetic field. The method 
we follow for impos ing random footpoint motion is de- 
scribed in detail in Longcope fl993) . We use a corre- 
lation time of Tc = lU. Only a band of Fourier modes 
with small k is driven so that the boundary flow is of the 
characteristic length scale of the order of the perpendic- 
ular simulation box size. The amplitude of the driving is 
small enough such that the root-mean-square boundary 
flow (urms ~ 0.075) is small compared with the Alfven 
speed along the large-scale magnetic field. 

The footpoint motion is assumed to take place on a 
timescale much longer than the characteristic time for 
Alfven wave propagation between 2 = and z = L, 
so that the plasma can be assumed to be in quasi-static 
equilibrium nearly everywhere, if such equilibrium exists, 
during this random evolution. For a given equilibrium, 
a footpoint mapping can be defined by following field 
lines from one plate to the other. Since the plasma is 
assumed to obey the ideal MHD equations, the magnetic 
field lines are frozen in the plasma and cannot be bro- 
ken during the twisting process. Therefore, the footpoint 
mapping must be continuou s for smooth footpoint mo- 
tion. Parker (Parker 1972) claimed that if a sequence 



of random footpoint motion renders the mapping suffi- 
ciently complicated, there will be no smooth equilibrium 
for the plasma to relax to, and tangential discontinuities 
(or current sheets) of the magnetic field must develop. 
Although Parker's claim has stimulated considerable de- 
bate, (see Janse, Low,& Parker 2010), we have shown 
elsewhere that it is valid if the equilibrium becomes un- 
stable because there is only one smooth equilibrium for 



a given footpoint mapping (Ng & Bhattacharjee 



19981 



under RMHD and using periodic boundary conditions 
in transverse coordinates. Moreover, thin current layers 
generally do appear in numerical simulations with small 
but finite resistivity, after random boundary ftows have 
been applied for a period of time. 

Ordinarily, the problem of calculating time-dependent 
solutions of Equations Q and ^ in line-tied magnetic 
field geometry involves all three spatial coordinates and 
time. As a first step in NB08, we had made a strong 
assumption that in addition to the coordinate z along 
which the magnetic field is line-tied, the dynamics de- 
pends on only one transverse coordinate x (as well as 
time t) . The nonlinear terms (those that involve Poisson 
brackets) in RMHD equations ^ and ^ then become 
identically zero. 

We have developed computer simulation codes that in- 
tegrate Equations ([ij and ([2| numerically for arbitrary 
footpoint displacements in both 2D and 3D. We use spec- 
tral decomposition in x and y (using a standard two- 
thirds de-aliased pseudo-spectral method) and a leapfrog 
finite difference method in z on a staggered grid. The ac- 
curacy of the Fast Fourier Transform (F FT) routine use d 
has been tested extensively as shown in Ng et al. (2008). 
For the 2D version, we can use an implicit method for 
time-integration so that we can take larger time steps 
than is allowed by the Courant-Friedrich-Lewy (CFL) 
condition for numerical stability of explicit methods, un- 
like in the 3D version, in which a predictor-corrector 
time-stepping is used. 



In NB08, the 2D code was run for cases with different 
resistivities. It was found that if Tc <^ tji, {tji is the char- 
acteristic resistive diffusion timescale), measured average 
rates of ohmic and viscous dissipation became insensi- 
tive to resistivity. A simple scaling analysis showed that 
this behavior could be derived beginning from general 
considerations if one includes the effects of the random 
walk nature of photospheric footpoint motions. How- 
ever, because the simulations lacked instabilities or mag- 
netic reconnection, the growth of magnetic energy was 
unbounded, with average transverse magnetic field By 



scaling as 77 



-1/2 



NB08 went on to demonstrate by fur- 



ther analytical argument that even if the growth of the 
magnetic field was thwarted, as would be the case in 3D 
simulations, dissipation would remain independent of re- 
sistivity, regardless of both the specific saturation level 
of B± and of the mechanism causing the saturation. It 
is this conjecture we seek to examine. 

The prescribed ordering Tc ^ te {te is the character- 
istic timescale over which energy is built up before being 
i mpulsively dissipated) is not met for the 3D simulations 
of Longcope & Sudan (1994), who found both dissipa- 
tion and By to scale as rf^^ . For our numerical experi- 
ment, we extend the range they examined by an order of 
magnitude in either direction, establishing an adequate 
domain to assess the stated hypothesis. We have also 
a dopted their foot point braiding algorithm (as described 
Longcope||1993 ) , allowing a direct comparison of scal- 



ing results in the range they examined. Before proceed 
ing with the details of our simulations, it is worth noting 
that NB08 demonstrated analytically as well as numer- 
ically that in the presence of steady footpoint motion, 
i.e., when (j)[z — 0) and (j){z = L) are time- independent, 
the heating rate is inversely proportional to ry, which is 
a strong and physically unrealistic dependence. For this 
reason, we do not pursue steady boundary flow, used in 
some other studies of coronal heating (e.g. Rappazzo et 



,al. „2008^ . 

3. NUMERICAL RESULTS: STATISTICAL STEADY STATE 

We have performed a series of simulations using our 
3D RMHD code described in Section [2] using a range 
of rj spanning two orders of magnitude to study scaling 
laws. Extending the range in rj by an order of magni- 
tude b eyond what was studied by [Longcope h Sudan 
( 1994 1 poses a significant challenge. As the dissipation 
coefficients (77 and v) get smaller, higher resolution has 
to be used to resolve smaller scales. To run the 3D 
code in high resolution, simulations are performed on 
parallel computers using MPI (Message Passing Inter- 
face). In order to run some cases for even longer time, 
we have also modified the code and run it on machines 
with CPUs (Graphics Processing Units) using Nvidia's 
Compute Unified Device Architecture (CUDA). 

The range of 77 has been extended to lower values (with 
Tc = 10 ^ Tr) for abou t an order of magnitude a s com- 
pared with the study in Longcope & Sudan ( 1994 ) which 
stopped at 77 ~ 10~^. This extension of course requires 
significant increase in resolution, with our highest reso- 
lution case at 512^ x 64 so far (runs with even higher 
resolution, such as RO in Table [l| have not been run for 
long enoug h time for good statistics), as compared with 
48^ X 10 inlLongcope & Sudani (|l994k. The main diffi- 
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culty in performing these simulations is the requirement 
to run up to hundreds or even thousands of Alfven times 
in order to obtain good statistics of the average quan- 
tities under the driving of random boundary flow. The 
basic parameters and results in this scaling analysis are 
summarized in Table [T] 

It is crucial to the scaling study that we obtain good 
statistics averaging over time evolution in statistical 
steady state. As with previous long time integration 
studies of the Parker model, the runs are started with 
a uniform magnetic field along e^. After initial tran- 
sients, the system will evolve to a statistical steady state. 
As mentioned above, thin current layers are formed and 
dissipated repeatedly during this statistical steady state. 
Figure [T] shows 3D iso-surfaces of J at a time taken from 
the R5 run, when there is a larger number of current 
sheets. This process is repeated indefinitely as the ran- 
dom boundary flows keep twisting the magnetic fields. 
Energy of the system is dissipated impulsively. As Poynt- 
ing flux injection progressively braids the fields, energy 
is built up until an instability drives current sheet for- 
mation and reconnection, after which energy is released 
in a short time. This is a major characteristic of the 
statistical steady state. For the runs R5 (Blue) and R12 
(Red), specified in Table [l] Figure [2] shows the intermit- 
tent nature of various quantities in time (in the unit of 
the Alfven time ta, the time it takes for Alfven waves to 
travel the distance of L = 1 between the two boundary 
plates along z). 

Figure l2r(a) shows the total magnetic energy Ej^i = 
jB^d^x, as well as total kinetic energy Ek — J v^d^x, 
where the integration is over the 3D simulation box. Note 
that the magnetic energy does not include the contribu- 
tion from the component, which is constant in the 
RMHD model. Since the applied photospheric fiow is 
chosen to be small (less than one tenth) compared with 
the Alfven speed (with Tc = IOta), the magnetic field 
configuration maintains quasi-equilibrium for most of the 
time, except when strong current sheets form from time 
to time, thus inducing instabilities and strong dissipa- 
tion. Therefore, Em is usually much larger than Ek- 

Figure [2] (b) shows the maximum current density Jmax 
over the whole 3D volume. Jmax increases over an or- 
der of magnitude on average in R5 compared with R12, 
and also fluctuates in time over a much larger amplitude. 
Observing Figures [2ja) and[2|b), note that the ratio of 
the increase in Jmax is much larger than the ratios of the 
increases of both Em and Ek as 77 decreases. 

Figure [2] (c) shows the ohmic dissipation W,, = 
77 / J^d^x. Similarly, there is also energy dissipated by 
viscosity, at a rate given by = ^ / ^'^d^x (not shown). 
Due to the same reason that Ek is much smaller than 
Em, the viscous dissipation is much smaller than the 
ohmic dissipation, if we choose v = rj (Prandtl number 
equal to unity), which holds for most of our simulations. 
In this case, the total energy dissipation rate (heating 
rate) is dominated by ohmic dissipation. When we use 
values of v greater than ry, the viscous dissipation can 
become a more significant fraction of the ohmic dissipa- 
tion. From the plot of Wrj, we note that even though it 
shows large fluctuations in time, it fluctuates around a 
higher level for R5 than for R12 due to the smaller value 
of resistivity in the former. 



To give a better measure of the level of energy dissipa- 
tion, we can calculate the time averaged energy dissipa- 
tion rates, e.g., Wr^ = [J* W,jdt']/t, and similarly for W„. 

The total energy dissipation rate is then W = Wrj + W^, 
which is plotted in Figure [2[d). Our physical assumption 
here is that such averagea quantities will tend to satu- 
rated levels as t tends to infinity. In practice, since we 
can only simulate for finite time, such saturated levels are 
found at a time t ^ ^ ta when these time-averaged 
values do iiot fiuctuate too much. We do see from this 
plot that W tends to saturate after a time much larger 
than TA- 

Also plotted in Figure l2jd) is the time average Poynt- 
ing flux /, where / — J v • Bd^a;, integrated over the 
top and bottom boundary surfaces with v = Up, the 
random photospheric flow. Note that / is not positive 
definite due to the fact that it involves the dot product 
between the velocity and magnetic field vectors and thus 
can be eit_her positive or negative. However, the time 
averaged / is almost always positive due to two factors. 
First, due to ohmic and viscous dissipation of energy into 
heat, if the total energy of the system is at a statistically 
steady level, there must be energy input from the bound- 
ary to provide this dissipation loss. Secondly, even when 
there is not much energy dissipation during a certain pe- 
riod, magnetic energy Em generally increases, since the 
magnetic footpoints at the two boundaries connected to 
the same magnetic field line will move apart from each 
other in a random walk fashion due to random photo- 
spheric motion. Therefore, a typical field line will gener- 
ally be stretched by the separation of the footpoints and 
the magnetic energy will increase. This increase in the 
magnetic energvmust corne from the Poynting flux. We 
see from FigureWd) that / tends to saturate in the long- 
time limit at a level close to that of W . In principle, these 
two rates should be the same, since the time averaged 
total energy also tends to a constant level. Numerically 
there is a slight difference between the two. Convergence 
studies show that this is mainly due to inaccuracy from 
finite grid size, and the difference decreases when higher 
resolution is used, especially in the parallel direction. 

Another measure of the accuracy of our solutions is to 
test the energy balance equation. 



d{EM + Ek) 
dt 



= I-W, 



V 



(3) 



Figure [2] (e) shows / as a function of time in pink for 
the run R5, d{EM + EK)/dt (calculated by taking finite 
difference in time) in purple, —W,^ — Wi, in blue, and the 
difference between the right and left hand sides of Equa- 
tion ([3]) in green. We do see that the residual power due 
to numerical inaccuracy is generally small compared with 
other terms. While accuracy can be improved by running 
at higher resolutions, doing so would require much longer 
run times, as well as limit the highest Lundquist number 
that can be simulated. In the context of energy balance 
in our simulations, we remark that the energy dissipated 
due to ohmic or viscous terms is essentially converted 
into thermal energy. No energy term or transport equa- 
tions are included. This is perhaps the primary weakness 
of this model, as it prevents us from predicting tempera- 
ture and density profil es which can be directl y compared 
with observations (see Dahlburg et al. 2009). However, 



Solar Coronal Heating 



5 



the heating rate required to maintain observed coronal 
tempera tures can ind e ed be estimated as has been done 
in, 



e.g.; 



Priest et al. (2002 1. NB08 followed this prac- 



tice and found that the heating rate determined from 2D 
simulations is consistent with such estimation, if the en- 
ergy dissipation does turn into heat as assumed. Readers 
should compare similarities and differences be tween this 



treatment with those used in other studies ( Longcope 



1993 



Longcope fc Sud an 1994 ; [Rappazzo et a. 



Hendrix & van Hoven 1996; He ndrix et al.||1996[ 



gaard fc Nordhmd 1996 ; Rappazzo et al.[[2010| T 

Figure[2](f) shows B± as a function of time. Here B± is 
defined as a root-mean-square value of the magnetic field 
strength, and so is effectively the square root of Em per 
unit volume. Similar to other time-averaged quantities, 
B± also tends to saturate at long time, at values that 
are much more physically reasonable than those in the 
2D runs. We remark that these saturated levels in the 3D 
runs are already much more reasonable than those in the 
2D runs that were found to have a scaling of B± a ri~^^'^, 
which can be much larger than unity (the value of the 
constant B^ used in the simulations). Thus, inclusion 
of 3D effects is seen to reduce the magnitude of B± to 
values that are smaller than the magnitude of B^- 

These effects include the formation of thin current lay- 
ers, onset of instabilities, and subsequent reconnection 
and enhanced energy dissipation. All these effects are 
more prominent when B± is larger, effectively limiting 
the growth of B±. Thus, these 3D effects can self- 
regulate the level of B± that can be built up when sub- 
jected to the driving of the random footpoint motion. 

Due to the fact that we are injecting energy into the 
system through random photospheric footpoint motion, 
a natural question to ask is whether this would induce 
other random processes, such as a turbulent cascade of 
energy that contributes to the heating of the corona. In- 
deed, turbulence has been studied in various coronal loop 
heating m odels beginning with the early work o f [van 



Ballegooije n 
drix et aliil 



( 1986) and more recently by others 



len- 



| |1996, Dmitruk et al.[[T998[ [Rappazzo et al 



2008p . However, as mentioned in the above discussion. 



we are driving with slow boundary flows (less than 1/10 
of the Alfven speed) with ta, and thus the magnetic 
field configuration maintains quasi-equilibrium most of 
the time. Moreover, we apply random boundary flows, 
instead of constant motion as in, e.g., [Rappazzo et al.[ 
(2008), so that energy injection is much slower due to 
the fact that magnetic fleld lines are stretched in a ran- 
dom walk fashion rather than at a constant rate. As a 
result, we have Ek <C Em, which is not consistent with 
equipa rtition of energy in A lfven wave turbulence. Sim- 
ilar to Hendrix et al. ( 1996 ) , energy spectra in our sim- 
ulations are largely exponential (not shown here) during 
relatively quiescent periods, with little or no impulsive 
energy release, but become progressively shallow power 
laws during particularly intense current sheet disruption 
events, with possible excitation of more Alfven waves for 
a short duration. As with their study however, compu- 
tation grid resolutions enable us to resolve less than a 
decade of the inertial range of the energy cascade. 

While it seems turbulence plays just a minor role in 
our present analysis, whether it plays a crucial role in 
determining the speed of magnetic reconnection has at- 



tracted a number of recent investiga tions (e.g. iLazar 
ian & Vishniac 1999; Lourc iro et al.|20 07; Bhattacharjee 
et al.|2009j |Loureiro et al.|2009[ [Kowal et al.||2009| ). It is 



evident there is a surge of interest m numerical experi- 
ments concerning turbulent reconnection and that much 
has yet to be settled. It would be interesting to see if any 
insights can be gleaned from our own data. As mentioned 
above, the presence of turbulence seems to be intermit- 
tent in our simulations, mainly during intense impulsive 
current sheet disruption events. Of crucial importance 
here is how well resolved we can be, and therefore how 
well developed an inertial range we can identify. This 
will depend on how low the value of rj (and thus how 
high the Lundquist number S) we can simulate in 3D, 
as well as how important physical properties scale with 
r/ or S, which we turn our attention to now. 

Figure [3] shows some of the scaling results we have ob- 
tained so far. In Figure [3] (a) , the time-averaged ohmic 
dissipation rate Wrj (at the saturated level) for different r/ 
for the runs listed in Table[l]^is plotted in triangles, while 
the viscous dissipation rate Wi, are plotted in squares. As 
pointed out above, W^, ^ in general, and thus the 
total dissipation rate (heating rate) W = -t- (plot- 
ted in asterisks) is very close to W^, except in the large 
resistivity limit. The time-averaged Poynting flux / is 
also plotted in the same graph in circles. It should be 
the same value as W theoretically, and we do see that 
the differences between these two quantities are gener- 
ally small in our numerical results, indicating acceptable 
accuracy. 

From this plot, we see that W actually only changes 
within an order of magnitude, and levels off at both the 
large and small 77 limit. This has important implications 
for the coronal heating problem, since the Lundquist 
number (on the order of the inverse of the normalized 
r/ in our simulations) can be as high as 10^** in the so- 
lar corona. Therefore, the leveling off of W at the small 
77 limit is especially important, and is in fact predicted 
by NB08 based on 2D simulations and theoretical argu- 
ments. As mentioned above, this level of W was shown in 
NB08 to be independent of dissipation mechanism pro- 
vided that the correlation time Tc is small compared with 
the time to build up magnetic energy. It was also esti- 
mated that this level of heating rate can give the same 
order of magnitude required for reali stic corona l heating, 
following similar considerations as in Priest et al. (2002). 
However, the amount of magnetic energy built up in this 
process does depend on dissipation mechanism, and be- 
comes unphysically large in 2D simulations in the small 
r/ limit (with Bj_ scaling as 77"^/^). We will now show 
that this scaling becomes much weaker in 3D. 

Figurep^(b) shows the time-averaged B±^ (at the satu- 
rated level) for different rj. This is a measure of the mag- 
netic fleld (or magnetic energy) production in the statis- 
tical steady state due to the applied random photospheric 
motion. Unlike W , B± production changes over an or- 
der of magnitude from large to small 77. This is because 
in the highly resistive limit, magnetic fleld produced is 
quickly dissipated and can only reach a low magnitude, 
while the dissipation rate does not decrease that much. 
In the small resistivity limit, the increase of B± slows 
down signiflcantly. 

Due to the fact that we are doing 3D simulations, and 
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that we need to simulate for a long time to obtain good 
statistics, so far we have only been able to extend the 
value of rj to about an order of ma gnitude lower, as com - 
pared with similar studies in (Long cope fc Sudari||1994 ). 
Nevertheless, we can see already that below ij ~ 10"'^, 
th ere is a significant deviation from the scalings obtained 



Longcope & Sudan (1994), who showed by numerical 



computation and scaling analysis that both W and B± 
should scale with ri~^/^ in the small rj limit. We have 
added dotted lines in Figure [s] (a) and (b) showing the 
77"^/'^ scaling. We see that the portion of the data in a 
range close to 77 ~ 10"'^ is indeed consistent with an r]~]/^ 
scaling. However, as described above, both W and B± 
increase much slower with decreasing 77 for even smaller 
77. This result has important implications on the solar 
coronal heating problem, since the Lundquist number in 
the solar corona is very high and hence we most likely 
need to have a mechanism to provide coronal heating 
that is independent of the Lundquist number in order to 
obtain a physically reasonable heating rate. At the same 
time, the magnetic field energy production should not 
increase to an unreasonable level compared with obser- 
vations. In addition to this numerical evidence, we will 
provide our own scaling analysis to make s ense of these 
results , as we ll as compare with results in [Longcope fc| 
Sudani (p94|. 



4. SCALING ANALYSIS 

We have shown an initial confirmation of the hypoth- 
esis of NB08 but as an additional goal we would like to 
understand the exact mechanism giving rise to satura- 
tion. Their conjecture made clear that the insensitivity 
to 77 holds true no matter what the saturation mecha- 
nism is. In order to provide a more complete numerical 
confirmation of their conjecture, it becomes necessary to 
identify the possible physical mechanisms behind satura- 
tion. 

A na tural place to beg in w ould be to examine t he re- 
sults of Longcope ( 1993 1 and Longcope & Sudan ( 1994 1 
who derived scaling laws based on Sweet-Parker recon- 
nection theory and analyzed a range in 77 which we have 
covered in our own study. By looking at where their scal- 
ing behavior or where their assumptions might be failing 
in our own numerical results, we might gain some insight 
into the physics occurring at even lower rj. The reader 
is referred to these papers for a detailed review of their 
scaling arguments. Here we will only discuss their as- 
sumptions and results briefly. 

They assumed Sweet-Parker theory is valid in the 
sense that when looking at only the current sheet region 
forming between two coalescing islands (flux tubes), the 
reconnection can be treated as a steady process in re- 
sistive MHD which results in the classic Sweet-Parker 
scaling relating the width 6 and length A of a reconnect- 
ing current sheet: 



6/ A - S 



-1/2 



(4) 



where S± = B^w/i] is the perpendicular Lundquist num- 
ber, with w being the perpendicular length scale of the 
reconnecting islands and so w ^ VpTc with Vp being the 
root-mean-square value of the random photospheric flow 
velocity. They also observed that both the number of 
current sheets N in the simulation box and the length of 



the current sheets A are relatively insensitive to resistiv- 
ity. We follow these assumptions as a starting point for 
our discussion, although we recognize that some of them 
need to be re-examined more carefully in future studies. 

In particular, the Sweet-Parker reconnection theory 
should apply only to higher Lundquist number (smaller 
rj) cases, in which the energy dissipation is dominated by 
the reconnection process. Therefore, the scaling analysis 
presented here should not work for larger 77 (i.e. 77 > 0.01 
here), which is actually not within the focus of our studies 
here. When the energy dissipation is mainly from the 
Sweet-Parker current sheets, the dissipation rate can be 
estimated by 



W ~ 77iVAL 



Bl 



BjLLl 

TE 



(5) 



where we have used the estimation that th^ current den- 
sity of the current sheet is given by J ~ B± / S and that 
the volume of the simulation box is LL^ . The energy dis- 
sipation timescale te in Equation ^ can then be solved 
as 

TE - L\/N{7jwB^f'^ (6) 

where we have used the Sweet-Parker scaling in Equa- 
tion Q. In a statistical steady state, the energy dis- 
sipation by Equation ^ has to be replenished by the 
production of magnetic neld energy due to the footpoint 
motion within the same amount of time te- 



In the st udies of Longcope ( 1993 ) and Longcope & Su 



dan (19941, although random photospheric motion was 
used in the simulations, the effects due to such random 
flows were not taken into account in their scaling anal- 
ysis. This can be justified if is much larger than the 
energy dissipation time te- In this case, the magnetic 
field strength production is given by 



Si 



B, 



VpTE 



BzVp 
LN 



WT] 



1/3 



(7) 



where we have used Equation (rol) and solved for B±. 
Putting back Equation ([T]) into Equation ^ results in 



TE 



1/3 



(8) 



W 



(9) 



^N'^wBzVpHj ^ 
and so the energy dissipation rate becomes 

(LfB^V" 

after putting Eqs. ^ and ^ into Equation ([s]). Note 
that all three of these quantities, B±, te, and W scale 
with ?7~ ^/'^, and thus we hav e recovered scaling laws de - 
rived in 
althoug' 



Longcope (1993) and Longcope &: Sudan| ( T994| , 
1 we are using a slightly ditterent approach! 
We may now put reasonable numbers into Eqs. ([t]) to 
([9]) and compare with our simulation results. Our sim- 
ulations are set up to use L = L± — B^ = 1- The 
root-mean-square photospheric flow velocity is measured 
numerically to be Vp ~ 0.075, and thus w ~ VpTc = 0.75. 
The average number of current sheets N is more difficult 
to determine and is subject to some uncertainties. How- 
ever, we have done some analysis (not shown here) of our 
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simulations for different rj and found that ~ 7 numer- 
ically in the small 7] limit. This seems to be somewhat 
higher than expected from the number of reconnecting is- 
lands (flux tubes). However, it is actually quite common 
to see multiple current sheets in a simulation output, as 
shown in Figure [T] 

Based on these values, we have te OJl/rj^^^, and 
so Te ^ 7.1 for 77 — 10^^. At the same time we get 
B_i_ 0.53 from Equation Q, and W 0.04 from Equa- 
tion ([9| at the same rj. Both of these are close enough 
to the values found in Figure |3] (a) and (b), and so it 
is an indication that our parameters used in these esti- 
mations are consistent with simulations. Compared with 
the value of Tc = 10, we see that although te is still 
smaller than Tc, it is about the same order of magnitude 
and thus Equation ([7| is only marginally justified. For 
larger rj, te is smaller, e.g., te ~ 3.3 for rj — 0.01 and 
thus is much smaller than Tc so that the random effect is 
not as important. This qualitatively explains why we see 
from Figure [s] (a) and (b) that there is a range roughly 
around rj ~ 0.01 to 0.001 where both W and B± scale 
approximately as as indicated by the two dotted 

lines in the two plots. However, for smaller rj, te becomes 
larger, e.g., te ^ 15 for r/ = 10"'* (if we continue to use 
the approximation, which may not be strictly valid) , and 
thus larger than te which makes the effect of randomness 
important. This explains the deviation from the 77"*/'^ 
scaling for both W and B± for rj smaller than around 
10-3. 

Now, taking into account the effect of random bound- 
ary flow, which makes the footpoints move in a random 
walk fashion as argued in NB08, the estimate for mag- 
netic field production must be changed from Equation ([7| 
to 



Vp{T,TEy^' 



L 



L 



N^wr] 



1/5 



^ (10) 
where we have again used Equ ation ([6|) and solved for 

B±. Substituting Equation (10 1 back into Equation (|6| 
results in 



TE 



L 



N^Tc \wBzVp7] 



-I 1/5 



and so the energy dissipation rate becomes 
W 



(11) 



(12) 



after putting Eqs. (10) and (111 into Equation ([5|), and 
thus it is independent of rj. Note that Equation (|12 1 is 
exactly the same as found in NB08 for systems reggirmess 
of dissipation mechanism, and is estimated to give the 
same order of heating consistent with observations. 

Using the same values of L, ij^, Bz, Tc, Vp, w, and N, 
Equation (11) becomes te ~ 0.42/r/^/^, and thus te ~ 
6.7 for J] = 10""^, if we could apply this equation. This 
turns out to be very close to te ~ 7.1 estimated above 
using Equation (|8]), indicating that the transition point 
between these two regimes of scalings is around 77 ~ 10^^ 
in our simulations. For rj = lO""*, Equation (111 gives 
Te ~ 17, which is significantly larger than Tc, and so 



these scalings based on random walk of footpoints are 
justified. 



Based on this set of parameters. Equation (12 1 pre- 
dicts W ~ 0.056 (independent of 77), which is close to 
the asymptotic values found in Figure [3] (a) in the small 
77 limit. We do see from this plot that W indeed does not 
increase as fast when ij is below 10""^, and is consistent 
with a trend to a constant level in small 7/, although we 
still only have a limited range of ri that we can simu- 
late. At the same time. Equation (10) gives a value of 
B± ^ 0.97, which is somewhat larger than expected from 
Figure [3](b), although we do need to recognize that there 
are uncertainties in these scali ng estimates. 

A better test of Equation ( 10 1 would be the scaling 
with rj in the small 7/ limit. In Figure [3] (b), we have also 
plotted a dashed line indicating the scaling of 77--'^/^. We 
do see that this seems to be consistent with a portion of 
the data of B± below 77 ~_10-'^. However, we cannot rule 
out the possibility that B± is actually increasing slower 
than r;-*/^, possibly due to a modification of the Sweet- 
Parker reconnection scalings, e.g.. Equation ([4|. We will 
further discuss this possibility in the next section. 

5. DISCUSSION AND CONCLUSION 

In this paper, we have presented our latest results 
based on numerical simulations of a 3D RMHD so- 
lar corona heating model for a range of 77 (and thus 
Lundquist number) with random photospheric motion. 
These simulations were performed over a period of more 
than two years and numerical results have been verified 
carefully to eliminate possible errors. So far, we have 
been able to simulate cases with rj about an order of 
magnit ude smaller than thos e presented in similar stud - 
ies in (Longcope 19931 and (Longcope & Sudan 1994 1. 
While this extension seems rnodest, it actually requires 
much more computational effort due to the increase in 
resolution and running time required, as well as the de- 
crease of time-step for numerical stability. To be able to 
achieve that, we have been running our simulations in 
parallel computers, as well as using CPUs. 

Moreover, we have shown that the extension of this 
scaling study towards smaller 77 turns out to have very 
important physical consequences. Numericjilly, we have 
shown that the scaling laws (wit h W and B± scale with 



-l/3^ 



) found in Longcope ( 1993 1 and Longcope & Sudan 



( 1994 1 become invalid tor 77 smaller than what was used 
in their studies (around 77 ~ 10"'^). Both W and B± 
are now found to increase much more slowly for smaller 
77, with W possibly leveling off to an asymptotic value. 
We have presented our own scaling analysis to justify our 
numerical results. By following s imilar assumpt ions as in 
Longcope (1993) and jLongcope fc Sudan (1994), e.g., us- 
ing Sweet-Parker scalings, we have been able to recover 
their 77" scaling laws for a range of 77 larger than 10""^. 
We have demonstrated that the transition between scal- 
ing behaviors derives from the fact that the effects of 
random photospheric motion are not important in the 
larger 77 range where the energy dissipation time te is 
smaller than the correlation time of the random fiow. 
For rj smaller than around 10"'^, te becomes comparable 
or even larger than Tc- In this range, an analysis based on 
the random walk of photospheric footpoints motion will 
predict the insensitivity to 77 we observe, further substan- 
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tiating the results found in NB08 which were based on 2D 
simulations and more general theoretical considerations. 
This is important to the problem of coronal heating since 
this heating rate has been shown to be consistent with 
the requirements for coronal hcati_ng. 

We have also shown that now B± has a much weaker 
scaling with 77, i.e., rj^^^^ instead. This is much bet- 
ter than the 77"^/^ scaling in 2D simul ations, as well as 
wea ker than the 77"^/^ scaling fo und in Longcope ( 1993 1 
and Longcope & Sudan (1994). This means that this 
scaling will predict a more physically realistic level of 
magnetic field as compared with observations. However, 
due to the fact that the Lundquist number of inverse 
of 77) the solar corona can be very high (up to 10^^ - 
10^'^), even a 77"^/^ scaling would result in an unrealis- 
tically large magnetic field, despite a much weaker de- 
pendence. The reason behind this is the fact that the 
Sweet-Parker reconnection rate, which scales with 77^^/^ 
is too slow for high Lundquist numbers. 

One solution for this problem is the possibility of a 
higher rate of magnetic reconnection even under resistive 
MHD. This possibil ity has attracted a number of recent 
investigations (e.g., Lazarian fc Vishniac 111999 Loureiro 



et al. 2007; Bhattacharjee et al.|2009 ; Loureiro etal.|2009 
I^wal et al.|2009ft. Many of these studies fall within tE^ 



scope of turbulent reconnection. While there are some 
indications that B± found in our simulations might actu- 



ally scale weaker than even 7/~^/^, we still have not been 
able to simulate even smaller rj to confirm this defini- 
tively. Moreover, the effects due to turbulence are still 
too difficult to study using our current level of resolution. 
However, this question is important enough that we are 
trying different ways to extend our range of rj to even 
smaller values to study these effects. 

In summary, by simulating with 77 about an order of 
magnitude smaller than in previous studies, we have been 
able to find new physical effects due to random photo- 
spheric fiows and thus new scalings with Lundquist num- 
ber. In future work, with further improvements in our 
computational approach, we hope to report results with 
even smaller values of 77, and investigate the possibilities 
of another asymptotic range involving secondary insta- 
bilities and turbulent processes. 
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TABLE 1 
Summary of Numerical runs 



Run 


V 








w 


W 


Poyiitj 111^ 


^ / ' A 


csol Lit ion 


RO 


0.00015625 


0.00015625 


0.542 


3470 


0.0446 


0.0103 


0.0587 


305.545 


1024^x128 


Rl 


0.00015625 


0.00015625 


0.537 


3440 


0.0468 


0.0127 


0.0546 


487.252 


512^x64 


R2 


0.00015625 


0.00062500 


0.610 


3900 


0.0513 


0.0283 


0.0519 


245.546 


512^x32 


R3 


0.00015625 


0.00062500 


0.614 


3930 


0.0498 


0.0275 


0.0458 


77.0269 


512^x32 


R4 


0.00031250 


0.00031250 


0.492 


1570 


0.0433 


0.00792 


0.0491 


857.407 


512^x64 


R5 


0.00031250 


0.00031250 


0.503 


1610 


0.0452 


0.00941 


0.0478 


9321.75 


256^x32 


R6 


0.00031250 


0.00062500 


0.502 


1610 


0.0431 


0.0111 


0.0467 


2032.02 


256^x32 


R7 


0.00062500 


0.00062500 


0.449 


718 


0.0416 


0.00540 


0.0427 


19342.8 


128^x32 


R8 


0.00062500 


0.00062500 


0.448 


717 


0.0399 


0.00502 


0.0401 


820.339 


128^x32 


R9 


0.0012500 


0.0012500 


0.372 


298 


0.0370 


0.00332 


0.0385 


11668.2 


128^x32 


RIO 


0.0012500 


0.0012500 


0.371 


297 


0.0373 


0.00336 


0.0411 


706.141 


64^ X 16 


Rll 


0.0025000 


0.0025000 


0.279 


112 


0.0299 


0.00272 


0.0311 


1317.70 


64^ X 16 


R12 


0.0050000 


0.0050000 


0.183 


36.7 


0.0215 


0.00317 


0.0252 


2566.96 


64^ X 16 


R13 


0.0100000 


0.0100000 


0.103 


10.3 


0.0132 


0.00394 


0.0168 


5209.60 


64^ X 16 


R14 


0.020000 


0.020000 


0.0547 


2.73 


0.00822 


0.00511 


0.0123 


10245.4 


64^ X 16 


R15 


0.040000 


0.040000 


0.0307 


0.767 


0.00623 


0.00544 


0.0113 


10240.3 


32^x64 


R16 


0.080000 


0.080000 


0.0197 


0.246 


0.00550 


0.00612 


0.0105 


10240.5 


32^x64 




Fig. 1. — 3D iso-surfaces of J at a time taken from the R7 run. Both figures are from the same time sample for the R7 run. The left panel 
shows iso-surfaces at J = —60, —24, 24, 60 while the right panel shows iso-surfaces at J = —36, —12, 12, 36. Blue and green iso-surfaces are 
made semi-transparent for greater visibility. 
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Fig. 2. — Plots (a), (b), (c), (d), and (f) show time series of various quantities for runs R5 (Blue) and R12 (Red). In (a), green and 
orange correspond to Ej^ for R5 and R,12 respectively while blue and red show E^. For plot (d), solid lines show W and dotted lines show 
I. For run R5 plot (e) shows -W = -(W,, + W^) (Blue), I (Pink), d{EM + EK)/dt (Purple; 
hand sides of Equation [s] (Green). Parameters used for these runs can be found in Table^ 
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Fig. 3. — (a) Average energy dissipation rate for different values of rj. A is ohmic dissipation, □ is viscous dissipation, O is the total of 
the two, and Q is the footpoint Poynting flux, (b) Average perpendicular magnetic field strength for different values of rj. 



